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ABSTRACT 


Bolton fc Schlegell (120101) presented a promising deconvolution method to ex¬ 


tract one-dimensional (ID) spectra from a two-dimensional (2D) optical fiber 
spectral CCD (Charge-Coupled Device) image. The method could eliminate the 
PSF (Point-Spread Function) difference between hbers, extract spectra to the 
photo noise level as well as improve the resolution. But the method is limited 
by its huge computation requirement and thus can not be implemented in actual 
data reduction. In this paper, we develop a practical computation method to 
solve the computation problem. The new computation method can deconvolve 
a 2D hber spectral image of any size with actual PSFs, which may vary with 
positions. Our method does not require large memory and can extract a 4k x 
4k noise free CCD image with 250 hbers in 2 hours. To make our method more 
practical, we further consider the inhuence of noise, which is thought to be an 
intrinsic ill-posed problem in deconvolution algorithms. We modify our method 
with a Tikhonov regularization item to depress the method induced noise, we 
do a series of simulations to test how our method performs under more real sit¬ 
uations with poisson noise and extreme cross talk. Compared with the results 
of traditional extraction methods, i.e., the Aperture Extraction Method and the 
Prohle Fitting Method, our method has the least residual and influence by cross 
talk. For the noise-added image, the computation speed does not depend very 
much on hber distance, the SNR converges in 2 to 4 iterations, the computation 
times are about 3.5 hours for the extreme hber distance and about 2 hours for 
non-extreme cases. A better balance between the computation time and result 
precision could be achieved by setting the precision threshold similar to the noise 
level. We hnally apply our method to real LAMOST (Large sky Area Multi- 
Object hber Spectroscopic Telescope, a.k.a. Guo Shou Jing Telescope) data. We 









find that the ID spectrum extracted by our method has both higher SNR and 
resolution than the traditional methods, but there are still some suspicious weak 
features possibly caused by the method around the strong emission lines. As 
we have demonstrated, our deconvolution method has solved the computation 
problem and progressed in dealing with the noise influence. Multi-hber spectra 
extracted by our method will have higher resolution and signal to noise ratio thus 
will provide more accurate information (such as higher radial velocity and metal- 
licity measurement accuracy in stellar physics) to astronomers than traditional 
methods. 

Subject headings: Algorithm; general — deconvolution, optimization, spectrum, 
optical fiber; individual — (LAMOST) 
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Introduction 


To extract ID spectra from 2D CCD images of optical fiber spectrographs, there are 3 
extraction methods that we often use. 

The hrst method is the Aperture Extraction Method (AEM, hereafter). This method 
sums up counts in CCD pixels within a given aperture along a hber trace on both sides of 
its trace center. It is quick and easy but fails to deal with cross talk between hbers when 
the distance is close. Meanwhile, simple summation will give higher weights to the pixels 
with lower Signal to Noise Ratio (SNR) at the wing of the prohle, which will deteriorate 
the whole extraction. 


To solve this problem, iHornel fllQSGI) and iRobertsonI fjl986l) introduced independently 
a method named the optimal extraction method. This method gives higher weights 
to pixels with higher SNR. After that, many people kept improving this method, for 


exampl^ 


traces. 


Lebonteiller 


Mukai 


f 19901 


Piskunov fc Valenti 


et ah 


Marsh 


1201011 


(1198911 . 


a pplied the method to extract spe c tra wit h disto rted 


Verschueren fc Hensbergel fll990l) . 


HvnesI ( 2002 1 and 


200211 developed many techniques to extract multi-hber CCD images. 


These technics are successful in some cases, yet they can not work correctly when there is 
cross talk between hbers. Because nois e of cu r rent C CD is very low, usually several counts 


Hornd (119861) . the optimal extraction method is 


per pixel, according to equation 11 in 
formally the same as AEM. 

With the development of the optimal extraction method, the Prohle Fitting Method 
(PFM, here after) arises. This method tries to simultaneously ht CCD counts of diherent 
hbers along the spatial direction with bell-like prohle functions at each wavelength. The 
spectral huxes can be derived from htting the parameters of the prohle fun ction. Thus, one 


can solve the cross talk 


(I2nnnh . 


Kelz et al 


pr oblem by httin g the 


(I2nn6h . 


overlap between hbers, e.g.. 


Ballester et al. 


SanchezI (2006) and lSharp fc Birchalll ( 2010 1 
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These traditional extraction methods all extract flux row by row, which assumes the 
flux distribution at every wavelength in a 2D CCD image is only a func tion of the spatial 


direction and independent of the 2D PSF (see 


Bolton fc Schlegell (2010|)). That is, flux 


I{x,y) measured at pixel {x,y) can be written as: 


X Hy), 


this equation indicates that the profiles along y direction at different x are not correlated, so 
that we can use the same ID profile to extract fluxes at different wavelengths. Considering 
the actual optical aberration, this requirement is too rigorous (see Fig. [T]). 


The fourth method is the deconvolution method, which was introduced by 


Bolton fc Schlegell fl2010l) . The principle of this method is that if we take a 2D spectral 


image as the convolution of PSFs and ID spectra, when the PSFs at all positions are 
known, the ID spectra can be obtained by deconvolution. This method can eliminate the 
profile difference between fibers, caused by the instrument optic al aberration. How e ver, t he 


biggest drawback of the deconvolution method, as discussed by 


Bolton fc Schlegell (120101) . 


is that it needs a calibration matrix made from PSFs, which is too huge to be stored in 
memory of common computers, leaving alone computation. In this paper, we will give a 
new solution scheme to solve this problem. We also try to apply the deconvolution method 
to simulations with noise and And that the method is very sensitive to noise, that is, the 
solution is very unstable when noise is induced. In this paper, we also describe how to 
partly solve the noise problem by adding a Tikhonov item in our method. 


This paper is organized as following: In Section 2, we introduce how to construct the 
PSFs from the arc lamp lines. Then in Section 3, our objective functions and algorithm are 
discussed. Computation tests on simulations with noise and serious cross talk as well as 
real telescope data will be illustrated in Section 4. The last section is our conclusion and 
discussion. 
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2. Constructing PSFs 

Due to instrument aberration, PSFs at different positions on a CCD image are 
different. To fully sample the actual PSFs on the CCD chip, non-blended emission lines 
at all positions on the CCD is necessary. This condition is hard to be satished without an 
optical frequency comb. Instead, we can sample basic PSFs from the arc lamp emission 
lines which are sparsely distributed along the wavelength direction. Since the PSFs change 
smoothly with position, for positions where there is no emission line, the PSF can be 
linearly interpolated from those basic PSFs. 

We select single lines from an arc lamp image. These arc lines should have high SNR 
and the distances between them should be as short as possible. Weak arc lines can be 
improved with another longer exposure. Then we calculate the centroid {x, y) of each 
good arc line, where y is the coordinate in the spatial direction and x is in the wavelength 
direction. The PSF of each arc line is normalized and then htted with a uniform B-spline 
surface. Then a basic PSF is gotten. 

Assuming PSFQ{xo,yo) is the PSF with center (a;o,|/o); PSFi{xi,yi) and FSF 2 {x 2 ,y 2 ) 
are two closest basic PSF neighbors on different sides of FSFq at {xi,yi) and {x 2 ,y 2 ) 
respectively. To calculate FSFq from FSFi and FSF 2 , we firstly move FSFi from {xi,yi) 
to {xi,yo), and FSF 2 from ( 0 : 2 , 1 / 2 ) to {x 2 ,yo), then get 2 new PSFs: FSF{ and FSF^. If 
FSF[i,j] denotes the PSF in the ith row and the jth column on the CCD image, then we 
can get FSFq by the formula: 


PSFQ[t,j] 


X2 -Xq 
X2 - Xi 


PSF[[t,j] + 


X2 - Xi 


( 1 ) 


As for the PSFs close to the image edge, we can extrapolate or replicate the PSF of 
the closet arc line. For the purpose of this paper, it’s not necessary to use complicated 
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algorithm, so the later is adopted. To save compuation memory, only the PSFs of arc 
lines are stored as basic PSFs during the program running, other PSFs are calculated by 
equation [T] whenever necessary. 

In some instruments, distances of neighboring hbers may be so close that PSFs can 
not be measured simultaneously in a single arc exposure due to £ber-to-£ber cross talk. 
This could be solved by masking all the other hbers on the focal plane while taking an arc 
exposure of a single £ber or hbers separated by large distances on the CCD chip at the 
same time. Combining PSFs from diherent exposures will give a full sample of PSFs on the 
CCD chip. However, the algorithm should not be applied if masking neighboring hbers is 
not possible for instruments with signihcant hber cross talk. 


3. The Objective Function and The Deconvolution Algorithm 
3.1. The Objective Fnnction 

Convolved with prohles of hber and instrument, a monochromatic light beam entering 
a spectrograph will end up with a bell shape PSF rather than a 6 function on the CCD. A 
2D hber spectrum can be treated as a lineup of these single PSFs at diherent wavelengths 
with diherent intensities: 


/ + 00 

C{x)PSFx,y{X-x,y)dx, (2) 

■oo 

where S{X,y) is the count recorded by the CCD at wavelength A and spatial coordinate y, 
and C{x) is the ID spectrum we want to derive. 

Considering cross talk between hbers in a multi-hber spectral image, the hnal count 
F{X,y) at position {X,y) can be written as: 


(3) 
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Nf 

a) = s.(A, </) + £,.„ 

i=l 

where Si{\,y) is the count of the ith spectrum at wavelength A and spatial coordination y, 
Nf is the number of hbers, and ex,y is noise. 


If we have an image of N pixels in the dispersion direction and M pixels in the spatial 
direction, then for an actual computer calculation, equation |3] can be written in a discrete 
form: 


Nf N 

F[/c, m] = jPSFij[k,m] + e[k,m], (4) 

i=i j=i 

where F[k,m\ is the flux at the fcth row and mth column of the CCD, Cij is the flux 
recorded at the zth hber and the jth row and e[k,m\ is the noise at the fcth row and mth 
column. 

From equation 01 we can see that it is a least-quare problem. We rewrite equation 0] in 
matrix form: 


F = AC+ e, (5) 

where A is the calibration matrix with FSFi^j[k, m] as its element at the {k* M + m)th row 
and the {i* N + j)th column, and C is a vector with qj as its (i * + j)th element. For 

convenience of description, we always assume that noise £ is independent and identically 
distributed in this paper, other cases can be easily changed into this case by a linear 
transformation. 


For gaussian noise, the objective function of equation 0] is: 
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N M Nf N 

[fc, m])\ 

k=l 111=1 i=l j=l 

For poisson noise, the objective function of equation [His: 


( 6 ) 


. ^ ^ {F% m] - ES Y.U m]f 


k=l m=l 


F[/c, m] 


We rewrite function [6] in matrix form: 


(7) 


- Acy {F - AC), 

C 


( 8 ) 


and its solution is 


(7= (V A)-^A^F. 


(9) 


If we denote W as a diagonal matrix with -p^:^ as its {k* M + m)th diagonal element, 
then we can rewrite function [7] as 


and its solution is 


min(F - AC)^W{F - AC), 
c 


C = (V WA)-^A^ WF. 


( 10 ) 


( 11 ) 


In the case of LAMOST f Cui et al 


(j2012|)), i.e., 250 spectra on a 4kx 4k image, if we 


directly solve equation [9] or [TT] as in the paper of 


Bolton fc Schlegell f 2010l) . the size of the 


calibration matrix A will be {Ak x 250) x {Ak x 4fc) = 1.6 x 10^^, so a double precision will 
require memory about 128TBytes. Even though matrix A is a sparse matrix, it can not 
be easily stored or computed. Therefore, we can not solve this problem directly by the 
least-square method. However, from function [6] and [7] we can see, they are both optimization 
problems. Therefore, we try to solve it by an alternative optimization method. 
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3.2. Algorithm 

For concision of expression, we will only discuss function [HI but the following procedure 
is also suitable for function [TOl 

As discussed above, direct computation of C from equation |9] will require unacceptable 
hardware resources. But if we begin with an initial guess of Cq, then change Cq in a way 
that could reduce the value of function [H] gradually, the final solution C will be found 
when the function converges to its global minimum. To further reduce the computation 
requirement, every time we only change a small block of C to minimize the value of function 
|8] while leaving other elements hxed. Then the next block of C will be changed to get 
another local minimum. At last, the global solution will be approached when the value of 
function |8] stops to decrease or the decrease is less than a given precision. 

To begin the iteration process, we first chose an initial guess, Cq, and a small piece 
of element block in Co, for example, the hrst k elements {k is a small integer). All 
the elements in this block will be changed to minimize objective function jH] Next, the 
\k/2 \ +lth,..., \k/2 \ +/cth element of the new Co will be selected as the new block, elements 
inside the block will be modihed to minimize the objective function while keeping other 
elements unchanged. Each time the block is shifted forward by [A;/2]. When all elements 
in Co have been updated, the new Co will be taken as the input in the next iteration. 

Since all PSFs in calibration matrix A are known, objective function [H] has only one 
global minimum, which is equation |9l If we keep repeating the above process, the objective 
function value will keep decreasing and achieve its global minimum point at last. 


R:=F- ACo 


We denote 


-11 - 

and 

res := R. 

For a given Co and a block made np of elements from the (s + l)th to the (s + A;)th, 
let 6C := (0,0, Xs+i,Xs+k, 0,0)'^, := (x^+i, x^+kV and define the (s + l)th, 

(s + k) th columns in A as submatrix A, then A6C = AXg. If we change Cq by (5C, the 
objective function becomes 

(F - kl(C'o + 5C)f{F - A{Co + 

={R- A5Cf{R- A5C) 

( 12 ) 

={R-AX,Y{R-AX,) 

=res - 2 R^Ax s + XjA'^AXs. 

If we let Xg = {A^A)~^A^R, where the objective function achieves its local minimum 
in the block Xg, then the value of the objective function becomes 

res — {R^A)(A'^A)~^{R^A)'^ < res. 

The second term in the left side of the formula is a non-negative quadratic term, this means 
the objective function will decrease when we change C locally by 6C. If we keep iterating 
in deferent blocks, the objective function will finally converge to its global minimum. 

The algorithm can be described as following: 

Step 0: Given a value e > 0, the size of block: k, and the size of variable G: N, 

0 < /c -C iV; 

Step 1: Let Co := 0, s := 0, ri := 0. For function [HI let R := F and res := R^R; 

Step 2: If s + k > N, s := N, goto setp 3; else let A be the submatrix consisting 
of the s + 1th ,...,s + fcth columns in A. We denote Xg := [AAA)~^AAR. Let 
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ri ri + 2R^AXg — X'^A'^AXg ; (c^+i,..., Cg+k}'^ (c^+i,..., Cs+fc)^ + Xg] R :— R — AXg 

and s := s + \k/2~\ + 1; goto step 2; 

Step 3: If s — /c < 0, s := 0, then go to step 4; else let A be the submatrix 
consisting of the s — k + 1th ,sth columns in A. We denote Xg := {A^A)~^A^R. Let 
rI -.= Ti + 2R^AXg — Xj'A'^AXg ; (c^-fc+i, Cg)'^ := {Cg-k+i-, ■■■■, + Xg] R := R — AXg 

and s := s — \k/2^ — 1; 

Step 4: res := res — ri] If ri < e, stop; else let ri = 0, goto step 2; 

We call this algorithm the Direct Deconvolution Algorithm, or DDA, hereafter. 

Since k is small, AXg in the algorithm can be sparsely stored and calculated. In 
principle, the algorithm can deal with an image of any size with the PSFs dependent on 
positions. 

In addition, for a multi-processor computer, if we use one processor to calculate a 
single hber at a time, this algorithm can be easily transferred into a parallel algorithm. 
For example, the size of our constructed LAMOST PSFs is 15 pixels x 15 pixels, if the 
distances between hbers are all greater than 8 pixels, then the odd hbers do not affect 
each other, neither do the even hbers. Then the following computation scheme could be 
applied: we hrst use processors to compute the odd hbers in parallel, after all the odd 
hbers have been computed, the even hbers will be computed in parallel. The above parallel 
computation process will be iterated until the solution is found. In general, if the minimal 
distance between the hbers is D, we can divide all the hbers into groups. We denote 

— G, and i = 1,2,..., 250 as the sequence numbers of the 250 hbers. The ^fth group 
consists of hbers with i mod G = g, where g = 0,1, 2,..., G — 1. Thus, hbers in the same 
group have no cross talk, our parallel algorithm can be applied in every group. 
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4. Computation tests 

Throughout this paper, the computer to test our algorithm is a MacBook Pro with OS 
Windows XP, Intel (R) Core (TM) i7, 2.67GHz CPU and 2G memory. The programming 
language is VG++ 6.0. 


4.1. Applied to a simulated 2D image based on LAMOST data 

To test our algorithm, we have constructed a simulation 2D image based on LAMOST 
data. 33 emission lines per hber from a LAMOST arc image are selected as basic PSPs 
to construct a 4A x 4iC GGD image with 250 fibers. 250 input ID spectra are extracted 
from a LAMOST target image by AEM while the spectral traces are derived from the 
corresponding flat field image. The size of each PSF is 15 pixels x 15 pixels. Each input 
ID spectrum (see Figj2]) is convolved with these PSFs along each spectral trace to generate 
a 2D spectrum for one fiber, then a 2D image with 250 hbers. 

In the computation process, the memory occupied is less than 300M Bytes, most of 
which is used for storing the original and residual images. To avoid singular matrix during 
computation when the viable block is too big, we set the size of the viable block to 20. 

Two tests are run following the steps described in section 13.21 the computation 
precision es are set to 10“^ and 10“^, respectively. After running the program for 7,074 and 
17,989 seconds, the results can be seen in FigJHl The top panel of Figj3] shows the input 
simulation image, the lower 2 panels show the residuals after the spectra are extracted with 
e = 10“^ and 10““^ respectively. For comparison, we put the 2D residual of the PFM in the 
fourth panel and the residual of the AEM in the bottom panel. From the corresponding 
color bars we can see that the residuals are extraordinarily small for our deconvolution 
method, and there are no obvious residuals around the emission lines in the 2nd and the 
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3rd panels. In contrast, the traditional PFM and AEM leave abont 5~7% residuals around 
strong emission lines. FigJU shows the absolute values of the differences between the input 
and the extracted spectrum. The upper 2 panels show the residuals of our deconvolution 
method with different e. Only one of the 250 spectera is shown to illustrate our results. We 
have checked all the other spectra, they are all at the similar residual level to the spectrum 
shown in FigJU As we can see, the residual is very small with about 2 hours computation 
when we set e to 10“^, and they are almost negligible but with a longer computation time 
when e is set to 10“'^. We also plot the residuals by the PFM and AEM in the lower two 
panels. Compared with the much larger residuals of the PFM and AEM, our method shows 
great advancement. 

As can be seen from the top panel of FigJU fhe absolute values of the residuals at the 
beginning part of the spectrum are relatively larger than the residuals after the hrst 500 
pixels. This is because that the calculation begins from one end of each spectrum in order 
of pixels. As described in 13.21 in each process of calculation, there is a block of \k/2\ pixels 
overlapping with the previous block. Results of these \k/2\ pixels in the previous block 
will be the input of the corresponding pixels in the current block. This guarantees that the 
solutions of conjunct blocks are related. Thus residuals in the later pixels will be smaller 
since their initial input had been processed closer to the hnal solution than the beginning 
part. The residual of the beginning part will be reduced if the program iterates long enough 
time, as shown in the second panel of FigJU 


4.2. Applied to simulations with cross talk 

In the above simulations, the distances between hbers were set to the same as 
LAMOST, about 15-16 pixels. The size of each PSF is 15 pixels x 15 pixels and the 
FWHM (Full Width at Half Maximum) is about 6 pixels, so £ber-to-£ber cross talk in the 
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above simulations is very small. To test our algorithm on images with serious cross talk, we 
construct 3 images using the same method discussed above but with the distances between 
fibers shrunk to about 10, 8 and 6 pixels, respectively (see FigJS]). We put 250 fibers in each 
image, the size in wavelength direction is kept as 4K pixels, so the image sizes are not the 
same for different fiber distances. Since the FWHM of LAMOST PSF is about 6 pixels, we 
do not try to simulate distances less than that. To study the influence of bright fibers on 
their neighbors, every 3 spectra, there is a spectrum with its flux 100 times stronger than 
its neighbors. 

The program runs about 2 hours by setting e to about 10“^, 2D and ID residuals are 
shown in Figj5]and FigJO] respectively. Only two of the spectra (a bright fiber and its faint 
neighbor) are shown in Figj6]to represent the common situation. From these two figures we 
can see: Firstly, both ID and 2D residuals are small and acceptable. The residuals increase 
when fiber distance decreases, but not in a linear tendency. In fact, the residual grows 
much faster as the distance decreases to the extreme case, 6 pixels. Secondly, residuals 
of the bright fiber and its faint neighbor are on the similar level and the influence of the 
bright fiber on its neighbor is negligible in our method. Thirdly, as discussed in Section 
14.11 residuals at the beginning parts of the ID spectra (Figj6]) are larger than in the rest 
parts. They could be reduced by setting computation precision e to a smaller value but 
with longer execution time. As shown in FigJTJ residuals at the beginning part are greatly 
reduced after the program runs about 5 hours when e is set to about 10“^. 


4.3. Comparison with Bolton & Schlegel’s method 


To compare with the method of Bolton & Schlegel, a piece of ID spectrum, with only 
1,000 pixels, is convolved with LA MOST PSFs to generate a 2D image of one fiber. Then 


we use the DDA and the method in 


Bolton fc Schlegell (120101) to deconvolute the 2D image 
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respectively. Results a re shown in FigJHl we can see that our method can achieve higher 


accuracy than that of 


Bolton &: Schlegell (120101) . When calculating an inverse matrix, 


calculation error is accumulated, so a larger matrix leads to a bigger calculation error. The 
calibration matrix (see equation [9] and IT^ in our method is much smaller than Bolton & 
Schlegel’s, so the calculation accuracy of its inverse can be higher than theirs. Meanwhile, 
calculation of large calibration matrix is very time consuming. It tooks about 123 seconds 
for Bolton & Schlegel’s method to run, while only 3 seconds for the DDA. Therefore, the 
DDA has better performance in both calculation accuracy and computation time. 


4.4. Applied to simulations with poisson noise 


In the above simulations, noise is not considered, yet in actual observations, noise is 
inevitable. Deconvolution is an ill-posed problem, its result is very sensitive to noise. From 
the upper 2 panels of FigUDlwe can see that if we use the DDA to direc tly extract an imag e 


with noise, th e 


Puetter et ah 


spec 


r um will be domina t ed by noise. Many people (e.g. 


Starck et ah 


(120051) , iDavies fc Kasperl (120121) a nd references the rein) 


( 120021 ) . 


to minimize the noise influence. Here we will use 


TikhonovI ( 19631 1 and 


rave discussed how 


Tikhonov fc Arsenin 


(119771) regularization methods to control the noise. Introducing the first derivative as 


regularization item in the Tikhonov method, function |6] could be rewritten as: 

N M Nf M Nf AT-i 

{F[k, m] — EE EE 

k=l m=l i=l j=l i=l j=l 


(13) 


where aij > 0 is a weight to adjust the regularization item (q^+i — Q^)^. We can rewrite 
equation [12] in matrix form: 


min(F - ACf{F - AC) + (FC)^a(FC), 
c 


( 14 ) 
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where T is the Tikhonov matrix consisting of 250 submatrix blocks on its diagonal. 

In each submatrix, the diagonal element {k,k),{k = 1,2,...,N — 1) is -1, the element 
{k, k + l),{k = 1,2,..., N — 1) is 1, and all the other elements are 0. a is a diagonal 
matrix consisting of 250 submatrix blocks, in the ith submatrix block, the diagonal element 
= 1,2, 1) is aij. 

For an initial guess Cq, we define the same variables: 6C, Xg, A, R and res as in 
section [3^ define the s + 1 th, ..., s + /c th columns in F as submatrix B, then A5C = AXg, 
V5C = BXg. Substituting the variables above into objective function (TU we could write 
the function in a new form: 

(F - A{Co + 5C)Y{F - A{Co + 5C)) + (F(C'o + 5C)YaiV{Co + 5C)) 

= {R - AXsf{R - AXg) + (FFo + BXgfaiTCo + BXg) 

(15) 

=res - 2RAAXg + XjA'^AXg + (FC'o)'^a(FC'o) + 2C^T^aBXs + XjB^aBXg 
=res + (FC'o)^a(FC'o) - 2{R'^A - C^T^aB)Xg + {XjA'^AXg + X^B^aBXg). 

If we let Xg = A + B"’"aB)~^{A^R — F^aFCo), where the objective function achieves 
its local minimum in the block Xg, and denote RES = res + (FC'o)^q;(FC'o), then the value 
of the objective function becomes: 

RES - {A^R - B^aTCofiA'^A + B^aB)-\A'^R - B^aTCo) < RES. (16) 

Similar to section 13.21 the second term in the left side of the function is a non-negative 
quadratic term, so the objective function will keep decreasing for every small block Xg in 
each iteration until the decrease is smaller than a given precision, the objective function will 
hnally converge close to its global minimum. Because F is a sparse matrix with non-zero 
elements all on the primary and the higher second diagonal, and a is also a sparse matrix 
with non-zero elements all on the diagonal, thus WaB and B^^qTCq can be sparsely stored 
and calculated. 
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During the computation process, denotes the decrement of the Tikhonov item in 
each iteration, ri denotes the decrement of res, other symbols are the same as dehned 
above. We then change the DDA slightly, and rewrite it as following: 

Step 0: Given a value e > 0, the size of block: k, and the size of variable C\ N, 

0 < A: < iV; 

Step 1: Let Cq := 0, s := 0. For function [TTl let R := F, res := R^R, F^ := 0, ri := 0, 
RES := res; 

Step 2: U s + k > N, s := N, goto setp 3; else let A be the submatrix consisting 
of the (s + l)th,..., (s + A;)th columns in A, and S be a submatrix consisting of the 
(s + l)th, ..., (s + k)th columns in F. We denote Xg = {A'^A + aB)~^{AAR — B'^aVCo)- 

Let n := ri + 2R^AXg - XjA'^AXg ; F^ := F^ - {2C^T^aBXs + XjB^aBXg); 

(cs+i,..., Cg+kV ■= (g+ 1 , Cg+kV + Xg; R := R- AXg and s := s + \k/2] + 1, goto step 2; 

Step 3: If s — fc < 0, s := 0, then go to step 4; else let A be the submatrix 

consisting of the {s — k + l)th ,..., sth columns in A, B he a submatrix consisting of the 
{s — k + l)th ,..., sth columns in F. We denote Xg = (A'^A + B^aB)~^{A^R — B^aTCo). 
Let n := ri + 2R^AXg - XjA^AXg; F^ := F^ - {2C^T^aBXg + XjB^aBXg); 

(^s—fc+i; •••: ^s) • (^s—fc+ 1 : •••; ^s) T R . R AXg and s . s [~A:/2~| 1, 

Step 4: res := res — ri and RES := RES — ri — F^; If ri + F^ < e, stop; else let 

ri = F(i = 0, goto step 2. 

We call this algorithm the Tikhonov Deconvolution Algorithm, or TDA, here after. 

To test the TDA, we generate a simulation image similar to that in Section 14.11 except 
that poisson noise is added to each pixel on the 2D image. Then the ID spectra are 
extracted by the DDA and the TDA. In the TDA, all the diagonal elements of a are set 
to 0.001. Results are shown in Figj9]and FigJTOl 2D residual images in Figj9]show little 
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difference between the DDA and the TDA, while compared with the DDA, the noise level in 
the ID TDA resnlt is greatly snppressed by the Tikhonov regularization, as show in FigJTUl 

We also extracted the noise-added 2D image by the AEM and the PFM. For 
comparison, absolute values of ID residuals of different methods are shown in Fig{TTl From 
the hgure, we can see, compared with the AEM and the PFM, residual of the TDA is 
smaller, so ID spectra extracted by the TDA are the most reliable. 

In the simulation image, the distances between LAMOST £bers(15-16 pixels) are the 
same as the size of the input PSF (15 pixels), then uncertainty caused by the hber cross talk 
is quite small, so the noise in the TDA ID spectra should mostly come from the method in 
extracting the 2D image with poisson noise. But since the TDA could extract the spectrum 
with the similar resolution to the input spectrum, whereas the traditional methods could 
not, so the large residuals in the AEM and PFM are partly from the resolution difference 
between the input and the output spectra. To compare the noise level of different extraction 
methods under the same resolution, both the input spectrum and the spectrum extracted 
by the TDA are degenerated to the resolution of the AEM by convolving with a gaussian 
prohle. Then the degenerated input spectra are subtracted from the ID spectra extracted 
by different methods. In FigJT^ the absolute values of residuals of the AEM, the PFM and 
the TDA are plotted in green, blue and red, respectively. From the hgure we can see, the 
residual of the TDA is obviously smaller than those of the AEM and the PFM. Therefore, 
our Tikhonov regularization item could successfully reduce the noise inhuence. 


4.5. Applying the TDA to simulation images with poisson noise and cross talk 

To test our method under more stringent conditions, we add poisson noise to the 2D 
images used in subsection 14.21 on which the hbers are separated by about 6, 8 and 10 
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pixels, respectively. We extract the noise-added images by the TDA. The weights of the 
Tikhonov item, a, are set to 0.0003 and 0.03 for the bright and the faint hbers respectively, 
the calculation precision e is set to 10®. The ID and the 2D residuals are shown in FiglT^ 
and Figim respectively. For comparison, the residuals of the PFM are also plotted in these 
hgures. We also list the corresponding SNR of each spectrum in FigJTdl here, the SNR 
is dehned as the mean of the extracted spectrum divided by the mean square root of the 
difference between the extracted spectrum and the simulation input spectrum. From these 
hgures we can see: Firstly, the 2D residuals of the TDA are comparable to poisson noise, 
which means that our can preserve exactly all huxes recorded by the CCD. Secondly, both 
the ID and 2D residuals of the PFM are larger than those of the TDA, the TDA spectra 
show much higher SNR than the PFM. Thirdly, as shown in FigJTTl the ID spectra of faint 
hbers are severely polluted by their bright neighbors in all the PFM results, meanwhile, the 
corresponding results of the TDA show that the residuals of both faint and bright hbers are 
around zero. The residuals in our TDA method show no strong increase with the decreasing 
distance between hbers. From these simulation results we can see, the TDA is much more 
reliable than the traditional methods especially when there is cross talk. 


4.6. Computation time 

The computation times in the above simulations are 167m, 168m and 252m for hber 
distances of 10, 8 and 6 pixels, respectively. As may be noticed, the calculation precision 
e is set to 10® in section 14.51 which is much larger than 10“^ in the noiseless simulation 
in section 14.21 As dehned in section 13.21 and 14.41 computation precision e is the threshold 
of the decrease of the objective function, which is summed over all pixels, that is, for our 
simulation with 250 hbers and 4k pixels in wavelength direction, the total number of pixels 
is much more than 10® (the accurate number depends on the hber distance, or how many 
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pixels the 250 fibers occupy), so e = 10® means we have an average difference much less 
than 1 per pixel . For the noise free simulation in section 14.21 our purpose is to test if 
our DDA method could work correctly with fiber-to-fiber cross talk, e is set to a small 
value to see if the result precision could be arbitrarily high when the calculation time is 
long enough. For the noise added simulation in section 14.51 the calculation precision is 
noise dominated, so e should be comparable to the noise level rather than the unnecessarily 
small value in section 14.21 We have made another test with e set to l.Oe — 20, the results 
show that both the 2D and ID residuals are at a similar level to that of e = 10®, but the 
computation time is 4 times longer. Beside hardware, computation time should depend on 
the image quality (e.g. noise, fiber number and fiber distance, etc) and input calculation 
parameters such as e, usually it is a case by case problem. It is better that we could know 
how to adjust the parameters to balance the computation time against the result precision. 
To better understand how the result precision changes with different conditions, we set e to 
an arbitrarily low value, i.e., l.Oe — 20, then for each iteration, we output the computation 
time, the change of RES (see section 14.4p . the SNRs of bright and faint fibers for the 
noise-added images of fiber distance 15, 10, 8 and 6 pixels respectively in Tabled! Only the 
first 9 iterations are shown in the table since the SNRs stop to increase for more iteration. 
From the table we can see, the computation time of each iteration are almost the same 
for different fiber distances, this is easy to understand since all the images we simulated 
contain the same number of fibers and the same size(4096 pixels) in wavelength direction, 
the amount of calculation is proportional to the product of the fiber number and the size in 
the wavelength direction, so the computation times are almost the same. For fiber distances 
of 15, 10 and 8 pixels, the result precision of larger fiber distance converges faster than the 
smaller ones, after the precision reach photon noise level (~ 10®), they all converge at the 
similar speed. But the precision decreases much slower for the extreme case, 6 pixels. If 
we use e = 10® as the threshold for our computation, the 15, 10 and 8-pixels converge at 
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a similar time, while the 6-pixels will take about 50% more time to reach the threshold, 
we highlight the computation times and the precisions in the table where they reach 10® 
threshold. In the right part of the table, we list the SNRs of one bright hber and its faint 
neighbor, the SNR is the same as defined in section 14.51 We have checked the SNRs of other 
hbers, the SNRs of faint hbers in the first 2 iterations change from hber to hber, which 
rehect the fact that the hrst result depends on the initial guess and is not stable, then the 
SNRs of diherent hber distances increase quickly to a similar level after the hrst 2 iterations. 
The SNRs for other bright hbers are similar to those listed in the table. From the table we 
can see that the SNR converges to a hx value even though the calculation precision keeps 
decreasing after each iteration, for the hber distances of 15, 10, 8 and 6 pixels, the SNRs 
converge after 3, 3, 4 and 6 iterations, respectively. The last column in Table [1] shows the 
SNRs of faint hbers with distance 6 pixels, we can see that the SNR increases much slowly 
after the 4th iteration, the corresponding precision is 6.89e8, comparing the precision and 
the computation time to those of 15, 10 and 8 pixels in the table, we can see that a better 
e to balance computation time against precision may lie between 10^ and 10^. From those 
tests we can see that for an image with noise, we could set e to a value that the average 
precision in each pixel does not exceed noise, usually 1 per pixel should be a lower limit. 
Once the precision threshold is set, except for the extreme cross talk, our method converge 
after 2 to 3 iterations, so the computation time does not depend so much on fiber distance 
for most cases. 


4.7. Applying the TDA to real LAMOST data 

Now we apply our TDA to extract an actual LAMOST image. FigJTBl shows the ID 
spectra extracted by the AEM(black curve), the PFM(green curve) and the TDA(red 
curve) from the same LAMOST 2D image, respectively. Most of the features in the spectra 
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are sky emission lines. We also plot a composite spectrum (blue curve) by convolving a 
gaussian pr o file (a =1 pixel, ~ 0.7A for LAMOST ) with the flux of sky emission lines from 


Hanuschik 


(1200311 ■ Relative fluxes and positions of those sky emission lines are indicated 


by vertical blue lines in the bottom. 


From FigJT5]we can see: 


Firstly, some emission lines indiscriminate in the AEM and the PFM results are clearly 
resolved by our method, for example, the strong emission line in the black curve around 
A6830 is clearly separated into two peaks in the red curve, which are OH A6829.5 and a 
blend line of OH A6832.7 and A6834.4. Other lines like A6842.0, A6871.1, A6894, A7041.5 
and emission lines between 6975A and 7000A are all clearly resolved by our method, as 
can be seen by comparing to the blue composite spectrum. Compared to the composite 
spectrum, the TDA result has similar resolution(R~4200), which is 3 times higher than the 
AEM and the PFM results(R~1000). 

Secondly, some ambiguous bifurcate peaks in the AEM and the PFM such as A6912.6, 
A6900.8, A6969.9 and A7003.9 are explicitly detected as single peaks by our method, which 
indicates that our result has higher SNR than those of the AEM and the PFM. 


Thirdly, carefully comparing the red spectrum with the black, green and blue spectra, 
we hnd some weak suspicious emission lines around strong peaks, such as A6855, A6918 and 
A6943, we think these features are most possibly caused by our Tikhonov regularization 
method, although the features of the underlying object can not be completely ruled out. 
Current Tikhonov item we adopted may produce wavy features around strong peaks, the 
influence of noise could not be completely eliminated. We will leave how to further constrain 
the noise influence in the deconvolution method for a future work. 
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5. Conclusion and Discussion 


In this paper, we present a practical calculation scheme to extract ID spectra from 2D 
images. Because of instrument optical aberration, PSFs at different positions of a CCD 
image are mostly different. To sample PSFs, we use emission lines on the arc image to 
generate discrete basic PSFs, then we use B-spline surfaces to represent the smooth basic 
PSF contours. By interpolating these basic PSFs, PSFs at non-emission-line region can be 
calculated. During the calculation process, only the smooth basic PSFs are stored in the 
memory, other PSFs can be calculated when necessary, so the requirement of memory is 
reduced. 


obiec 


Due to its huge requirement of computation resources, it is im possible to solve th e 
ive function directly by the least-square method described in 


Bolton fc Schlegel 


(120101) . To solve this problem, we try to reduce the objective function in a small variable 


block in one calculation step, thus only a small amount of memory is needed. The objective 
function could converge to its global minimum by gradually decreasing the variable blocks 
to their local minimum. Thus, our calculation scheme can solve the problem with common 
calculation resources. 


Based on our calculation scheme, we apply the Direct Deconvolution Method to extract 
simulation 2D images without noise. Results show that our methods could extract the ID 
spectra in a precision that traditional methods could not achieve, even if fiber-to-fiber cross 
talk is significant. 

Deconvolution methods are sensitive to noise. To suppress the noise influence, we 
introduce a Tikhonov item into the objective function. Then we apply the TDA to 
simulations with poisson noise. Results show that the TDA can extract the most reliable 
ID spectra compared with the AEM and the PFM. It can correctly extract the spectrum 
with extreme cross talk, even for the faint fibers with fluxes 100 times lower than their 
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neighbors. 

We also apply the TDA to an actual LAMOST 2D image. Compared with the AEM 
and the PFM, the TDA can extract ID spectra with both higher SNR and resolution. 
Theoretically, our algorithm discussed in this paper can be applied to a 2D image of any 
size and PSFs with any shape, even there is serious cross talk between hbers. 

We further consider how to set precision parameter e to balance the computation time 
against the result precision, and hnd that setting the calculation precision similar to the 
noise level should be a reasonable choice. Because our method converges very fast, for 
images with not so extreme cross talk, the computation time does not depend so much on 
hber distance, usually in 2 to 3 iterations, the result will be around the peak SNR. For our 
simulation with 250 hbers and 4k pixels in wavelength direction, the computation times are 
about 2 hours for the images with hber distance greater than 6 pixels and 3.5 hours for 6 
pixels. Moreover, our algorithm could be easily transferred to a parallel algorithm, and the 
computation speed will be boosted. 

In summary, compared with the traditional PFM and AFM, deconvolution method is 
the most consistent with physical process that 2D spectra are recorded by the CCD, so it 
can extract the most accurate huxes from a 2D image, especially for multi-hber spectra. 
Beside these, spectra extracted by our deconvolution method have higher resolution and 
SNR, so it is the most promising extraction method. 


However, there are still some uncertainties on weak features around strong emission 
lines in the ID spectra extracted by the TDA. Laborious ehorts are still needed to solve 
the problem of noise inhuence. Besides the TDA, there are many other methods in the 
literature to attenuate the noise influence in signal processing, such as the linear regularized 
methods, the Baye sian methodology and the wavelet-based deconvolution methods (see 


Starck et al. 


(120021) 1. To study how to attenuate the noise influence on the deconvolution 





method will be the topic of our future work. 
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Table 1: Computation time and calculation precision in each iteration for images with different fiber distances 



Computation time(s) 

Precision (difference of RES, see sectiorl4.4n 

SNR(Brihgt Fiber) 

SNR(Faint Fiber) 

fiber 

distance 

Iteration 

15 

10 

8 

6 

15 

10 

8 

6 

15 

10 

8 

6 

15 

10 

8 

6 

1 

2487 

2467 

2520 

2472 

2.35el3 

2.35el3 

2.35el3 

2.30el3 

3.58 

3.55 

3.57 

3.62 

4.57 

1.38 

0.17 

0.04 

2 

5035 

5015 

5030 

5032 

2.95e9 

3.50e9 

3.58el0 

6.28ell 

15.07 

14.74 

15.04 

14.16 

4.58 

4.74 

4.19 

0.80 

3 

7550 

7507 

7557 

7535 

2.12e7 

2.13e7 

7.89e7 

2.01el0 

16.34 

16.18 

16.70 

16.06 

4.58 

4.71 

4.71 

4.16 

4 

10057 

9997 

10085 

10027 

2.16e5 

2.62e5 

3.65e5 

6.89e8 

16.27 

16.14 

16.68 

16.02 

4.58 

4.71 

4.70 

3.98 

5 

12555 

12490 

12602 

12575 

4.27e3 

4.37e3 

4.50e3 

2.39e7 

16.25 

16.13 

16.67 

16.00 

4.58 

4.71 

4.70 

4.27 

6 

15057 

14977 

15087 

15120 

84.9 

84.9 

87.2 

9.55e5 

16.25 

16.13 

16.67 

16.00 

4.58 

4.71 

4.70 

4.30 

7 

17552 

17470 

17600 

17662 

1.94 

1.93 

1.99 

4.15e4 

16.25 

16.13 

16.67 

16.00 

4.58 

4.71 

4.70 

4.31 

8 

20097 

20017 

20092 

20167 

4.90e-2 

4.89e-2 

5.05e-2 

1.86e3 

16.25 

16.13 

16.67 

16.00 

4.58 

4.71 

4.70 

4.31 

9 

22592 

22570 

22590 

22707 

1.35e-3 

1.35e-3 

1.40e-3 

7.75 

16.25 

16.13 

16.67 

16.00 

4.58 

4.71 

4.70 

4.31 
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ig. 1 .— The left panel is a piece of arc lamp fiber spectral image of LAMOST (jCui et al. 


201211 1. bright spots are emission lines. The green and the red lines indicate rows where we 


sample the ID prohles on the same arc emission line. The right panel illnstrates normalized 
prohles sampled on the row with the same colors in the left. The prohle difference is obvions. 
We can see that extraction method based on I{x) x I{y) can not work properly when the 
prohle is position dependent. 
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Fig. 2.— One of the input ID spectra used to generate the 2D simulation image. It is the 
red part of a LAMOST spectrum. As we can see, there are many strong OH sky emission 
lines. These lines will help us to evaluate how our deconvolution method works on reliability, 
SNR and resolution. 
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796.32 

637.06 
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318.53 

159.26 

0.00 



1.43E-012 
9.14E-013 
3.97E-013 
-1.19E-013 
-6.35E-013 
-1.15E-012 
-1.67E-012 
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4.18E-015 

-1.89E-013 

-3.83E-013 

-5.76E-013 




123.30 

88.11 
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17.73 

-17.45 

-52.64 

-87.83 




Fig. 3.— The top panel is a piece of simulation image, we apply different methods to this 
simulation image, the residuals are shown in the lower 4 panels. From the 2nd to the bottom 
panels are the residuals of the DDA algorithm with e = 10“^, the DDA with e = 10“*^, the 
PFM and the AEM, respectively. Different colors represent count level on the CCD image as 
indicated by the color bars on the right. Compared with the PFM and the AEM, residuals 
of the deconvolution method are extremely small, especially around the region of emission 
lines. 
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Fig. 4.— This figure shows the absolute values of the difference between the input spectrum 
and the spectrum extracted by different methods. From top to bottom, the methods are: 
the DDA with e = 10“^, the DDA with e = 10“^, the PFM and the AEM, respectively. 
Comparing with the traditional PFM and AEM methods, we can see that the residuals of 
deconvolution method are negligible. 
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1.2e-011 
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-8.3e-013 
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1.2e-008 

4.0e-010 

-1.1e-008 

-2.2e-008 

-3.4e-008 
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-0.068 

-0.14 

- 0.22 


Fig. 5.— From top to bottom, the fiber distances in each row are set to 10, 8 and 6 pixels, 
respectively. Panels on the left colnmn are pieces of simnlated 2D images. Every 3 spectra, 
there is a spectrnm 100 times stronger than its neighbors, only the bright fibers can be seen 
dne to high contrast. Panels on the right colnmn are the corresponding residnals after the 
algorithm rnns for abont 2 honrs. Connts are represented by colors as indicated by color 
bars. 







































35 


C 

o 

CJ 



1000 


2000 

pixe 


5000 


4000 


c 

o 

o 



pixe 


1 0 ^ 

10 ' 

10° 

^ 10^' 

8 10 -^ 

1 0^° 

1 0 ^^ 

1 0 ^° 

0 1000 2000 5000 4000 



Fig. 6.— Absolute values of ID residuals after the algorithm runs for about 2 hours (e ~ 
10“^). From top to bottom, panels illustrate absolute values of residuals of two ID spectra 
extracted from 2D images with hber distance about 10, 8 and 6 pixels, respectively. In each 
panel, the red spectrum is the residual of the bright hber and the black spectrum is its 
neighbor with hux 100 times weaker. 
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Fig. 7.— The same residuals as FigJHl but for the program running for about 5 hours 
(e ~ 10-^). 
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Fig. 8.— The top panel is the input ID spectrum used to generate 2D image; the middle 
panel is the difference between the input spectrum and the spectrum extracted by Bolton & 
Schlegel’s method; the bottom panel is the difference between the input spectrum and the 
spectrum extracted by the DDA. Comparing the lower two panels, the accuracy of our result 
is much higher than Bolton & Schlegel’s. 






38 




973.00 

810.83 

648.67 

486.50 

324.33 

162.17 

0.00 




80.60 

54.05 

27.50 

0.95 

- 25.60 

- 52.15 

- 78.70 




80.60 

54.05 

27.50 

0.95 

- 25.60 

- 52.15 

- 78.70 


Fig. 9.— The top panel is the simulation 2D image with poisson noise, the middle panel is 
the residual of the DDA and the bottom panel is the residual of the TDA. Pixel counts are 
represented by different colors, as indicated by color bars on the right. The mean of the DDA 
and the TDA residuals are -0.0028 and 0.0026 respectively, and the standard deviations of 
the DDA and the TDA residuals are 12.23 and 12.75 respctively, so the 2D residuals of the 
two methods are almost the same. 
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Fig. 10.— From top to bottom: the first panel shows the input spectrum (red curve) and the 
spectrum extracted by the DDA from the noise-added image shown in Figj9](blue curve); the 
second panel shows the difference between the two spectra in the hrst panel, from which we 
can see the influence of noise on the DDA is enormous; the third panel shows the spectrum 
extracted by the TDA (blue curve) and the same input spectrum as in the top panel (red 
curve); the last panel shows the difference between the two spectra in the third panel. 
Comparing the second panel with the fourth panel, noise introduced by the DDA is greatly 
suppressed by the Tikhonov regularization. To make a close view of the spectra, only a part 
of the full wavelength coverage is displayed here. 
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Fig. 11.— To compare the residual of the spectra extracted from the noise-added image by 
different extraction methods, we show the absolute value of the AEM residual (green curve) 
in the upper panel and the absolute value of the PFM residual in the lower panel (blue cure), 
the absolute value of the TDA residual are plotted in red in both panels. Compared with 
the AEM and the PFM, the residual of the TDA is the smallest. 
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Fig. 12.— The same residual plots as FiglTTl but to compare the noise level of different ex¬ 
traction methods in the same resolution, both the input spectrum and the spectrum extracted 
by the TDA are degenerated to the resolution of the AEM and the PFM by convolving with 
a gaussian prohle. Compared with the AEM and the PFM, the influenced of noise in the 
TDA is much less even if the resolution is the same. 
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Fig. 13.— From left to right, columns are 2D residuals of the PFM, 2D residuals of the TDA 
and the input 2D poisson noise,respectively. From top to bottom, the distances between 
hbers in different rows are 10, 8 and 6 pixels, respectively. Comparing the medium to the 
right column, the 2D TDA residuals are at the similar level of the poisson noise. 
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Fig. 14.— ID residuals from images in FigJTSl The left and right columns show ID residuals 
of one of the bright hbers and its faint neighbor respectively. From top to bottom, the 
distances between hbers are 10, 8 and 6 pixels, respectively. In each panel, the red spectrum 
is the residual of the TDA and the black is the residual of the PFM. We also calculate the 
average SNRs of the extracted spectra, which are marked with the corresponding colors in 
each panel. 
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Fig. 15.— Results of extracting an actual LAMOST image. The spectra extracted by the 
AEM, the PFM and the TDA are plotted in black, green and red, respectively. The blue 
spectrum is a composite spect rum by convolving a gaussian prohle with sky emission lines 


fluxes from 


Hanusc 


lik 


and positions from 


(120031) . We also plot the sky emission lines with relative fluxes 


Hanuschik 


(120031 ) in the bottom. Each spectrum is scaled and offset 


arbitrarily in y axis for clarity. 




























































